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The influence of dilution with non-magnetic impurities on the order in the ground state of the 
ON S = 1/2 Heisenberg antiferromagnet is investigated by Quantum Monte Carlo simulations. Data of 

On . the spin correlation functions and thermodynamic properties for system sizes up to L x L = 16 x 16 

and for impurity concentrations up to 5 — 37.5% are presented. In the low doping regime the 
^-j, correlation length shows similiar dependence on temperature as in the pure case. The staggered 

magnetization (N^(S)} in the ground state which is the order parameter of the present model is 
estimated by extrapolating data at low temperatures. Impurity concentration dependence of {N^(S)} 
is presented. (N^(S)} decreases monotonically with impurity concentration and becomes very small 
around S c ~ 0.35 which is much smaller than the classical percolation threshold. 

PACS numbers: 75.10.Jm, 75.30. Kz 
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In this article we show new results on the spin correlation functions and the long range order of the diluted 
Heisenberg antiferromagnet on the square lattice with S — 1/2. Originally this work was intended to understand the 
effects created upon doping of the high temperature superconductor materials.. It is now generally believed that the 
pure antiferromagnet with S = 1/2 on the square lattice has long range ordeiufi. But the value of the order parameter 
is much reduced against its classical value due to quantum fluctuations. The interesting point now is the influence of 
C | inhomogeneities on the competition between order and quantum fluctuations. 

. For our purpose it is necessary to study large lattices and these systems are too large to use exact diagonalization. 
Therefore we used the Quantum Monte Carlo method (QMC). Because we are working only with static impurities, 
the minus sign problem does not suffer the present QMC. And our results might be applicable to experiments with 
O ' materials like LaiCu\^^Mg^O^ and La2Cui-sZnsOjj- In these materials magnetic copper atoms are substituted by 
, non-magnetic zinc or magnesium atoms, so there are static non-magnetic site impurities, which are the same type of 
? 1 impurities we are considering in our simulations. 

r> | In the present paper we report the temperature and impurity concentration dependence of the spin correlation 
. functions and also the long range order in the ground state. 

We have studied the spin correlation functions for various concentrations in a wide temperature region including 
lower temperatures than the previous Manousakis' worka. We analyze the concentration dependence of the correlation 
length in the form 

£(T,5) = Q(«5)e^ . (1) 

The correlation data fit this behavior well. And the estimated spin stiffness ps(5) is monotonically decreasing with 
impurity concentration. 

In the classical analog of the model, the long range order exists as long as the lattice itself percolates. Namely, the 
long range order vanishes at an impurity concentration corresponding to the percolation threshold, 6 C = 0.41B. In the 
present paper, the impurity concentration dependence of the long range order is reported. We would like to point 
out that the amount of the long range order becomes very small for 6 ~ 0.35, which is much smaller than S c — 0.41. 
The concentration dependence of percolating lattice sites has a dependence (S — S c )" with j3 < 1 and therefore it has 
a convex shape. On the other hand, the obtained dependence of the long range order for the present model has a 
concave shape. Furthermore, it is suggested that the critical value of the concentration is also reduced in the quantum 
case. 
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Our method is described in chapter 0. The results for the spin-spin correlation functions are shown in chapter III 
and those for the long range order in chapter IV. The results are summarized and discussed in chapter |v|. 

II. METHOD 

The model including the inhomogeneities used here can be described by the Hamiltonian 

n = jJ2 £ i E A ■ &i > ( 2 ) 

<i,j> 

where Si = for holes otherwise e» = 1 (for spins), Si are the S — 1/2 spin operators (Si — \<?i), and the sum 
runs over nearest neighbour pairs i and j, J is the antiferromagnetic coupling constant. The positions of the holes 
were chosen randomly but fixed for one simulation, so they are static non-magnetic impurities. The concentration of 
impurities is defined by S, 

5 = l-^>/L 2 , (3) 

i 

where L is the linear system size of the square lattice and N = L x L is the number of lattice sites. 

For the simulatioiwtt£e. used the Quantum Monte Carlo method based on the Suzuki-Trotter formula and checker 
board decompositiorooEl. Periodic boundary conditions were applied to the x-, y- and Trotter directions. In addition 
to standard local flip types we used special flip types in the neighbourhood of the impurities, local loop flips around 
small clusters o£,holes or single holes and non straight global flips in the x-y plane. Details of these new flip types 
are described inH. These new flip types are necessary for the ergodicity of the algorithm. 

We study system sizes L x L with L = 4, 8, 12, 16 for impurity concentrations up to S = 37.5% corresponding 
to 6, 24, 36, 96 holes. We perform simulations at various temperatures. The lowest one is T = 0.05J. For this 
low temperature it is necessary to go up to 9 • 10 5 Monte Carlo steps and Trotter numbers m = 32. We perform 
simulations for each sample with 3 or 4 different Trotter numbers m and then extrapolate for m to infinity by using 
a linear function in (^r). An extrapolation up to quadratic order is used as a test for the error of the extrapolation. 
From these data the ground state properties are obtained by extrapolations. 

The data for L — 4 are available from exact diagonalization. From these data we estimate the low temperature 
properties by summing up several low lying levels with the canonical weight. They are in good agreement with the 
data obtained by QMC. Larger systems are simulated by QMC. For each system size and hole concentration we 
average the data over five different samples with hole positions chosen randomly. Between them we found only small 
variation. Averaging over more than five samples is too much computer time consuming. 

To study such a big parameter regime it is necessary to make an optimized and fully vectorized program. Locations 
of the new type flips for a given hole configuration are listed up and stored before the simulation. All transition 
probabilities for each fliptype, which are given by ratios of Boltzmann factors between possible local spin configurations, 
are also stored in tables in the beginning of the run. During the simulation the computer only has to compare the 
transition probability for a flip with a random number. If the flip is accepted, the configuration is replaced by the 
flipped configuration which is also stored in another table. 

After one sweep, namely after each flip type is tested over the whole lattice, physical quantities are calculated. The 
data of 10000 sweeps are combined to one bin. The samples for each bin seem to be statistically independent with 
Gaussian distribution, so they can be used to calculate statistical errors. These data are extrapolated for m to infinity 
as mentioned above and then averaged over different configurations with the same system size and hole concentration. 
Averaging and extrapolation is done using a weighting procedure, each data point is weighted by the reciprocal of the 
square of its error. Final error bars are estimated only by this procedure. The further treatment of the correlation 



data are described in chapter III and the sublattice magnetization in chapter IV 



III. CONCENTRATION DEPENDENCE OF THE SPIN CORRELATION FUNCTIONS 

One interesting property of the system is how the spin correlation function depends on the distance between two 
spins 

<r) = (-iy^J2{S?S? +r ), (4) 

i 
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where we only change the positions along the x- or ^-direction. 

In Fig. (]]]) the spin correlation functions for the system size L = 16 and temperature T = 0.1J is shown for 
different hole concentrations. The correlation functions decrease with distance, although the data are symmetric at 
r = 8 because of the periodic boundary conditions. They also monotonically decrease with hole concentration for each 
distance. In referencescHlju the enhancements of short range correlation function near impurities are reported. But 
in the normalization of Eq. (Q) the enhancements by nearest neighbour correlations near impurities are compensated 
by decreasement of the number of pairs of spins. 

In order to study the temperature dependence of the correlation length, we have to fit our data to the expected 
correlation function. There is no agreement in the literature about the right form of this function as can be seen in 
the discussion 00. Ding and Makivi<x3 used the form 

c(r) = ylr~ A e~f , (5) 

where they used both of A and A-as the fitting parameters and found A ~ 0.4. However spin wave calculations!^ 
and the-S|ehjsvinger boson methodEfl give A = 1, the renormalization group approacbo gives A = 0.5 and some other 
authorsEjIIja do not use the algebraic decaying part corresponding to A = 0. _ 

If we generalize the temperature dependence of the correlation length proposed by Chakravarty et. al.El to the case 
with impurities, we haveo 

£(T,5) = Q(£)e^ . (6) 



The prefactor C{ and the spin stiffness ps now depend on the hole concentration 5 and we want to test the validity 
of the form and specify the parameter regime where the form is valid. 

We tried the least square fits of the data of the correlation function to the form Eq. (^) with different values of A. 
In Fig. (0) the fitting curves using A = 1/2 are shown. The fits look excellent and the errors for the fitted correlation 
lengths are very small (if we exclude the nearest neighbour spin pairs for the fit). If we use A = 1 we get nearly the 
same results if the fitted correlation length is smaller than the system size. We find no significant difference for these 
two values of A. But if we put A = we find bigger deviations of the fitting functions from the data points, and the 
fitting errors are much enhanced. In our analysis, there are too small number of data points to use A also as a fitting 
parameter. If we would use the two parameters, A and £, as the fitting parameter, A varies largely to fit the shapes. 
In particular, for the cases where we found a long correlation length by using a fixed value of A, the two parameter 
fitting concludes an algebraic decay without exponential part, namely A 7^ and £ ~ 0. 

So we estimated the correlation lengths by fitting using A = 1/2. Their logarithmic values are plotted in Fig. (||) 
over 1/T for the different hole concentrations. Only the correlation lengths which are smaller than the system size 
are plotted. For these data the size dependence is small. Data appear in straight lines in this graph. Thus the 
temperature dependence is consistence with Eq. (fy and we can estimate the spin stiffness ps from the slopes. In this 
graph error bars are not drawn because some of them are overlapping and they are rather confusing for the eye. For 
temperatures higher than T — 0.5J the error bars are nearly like the point size. At T — 0.25J, the error bars for 
low concentrations, where the correlation length becomes comparable to the system size, increase up to 10%. In the 
estimation of ps, this fact is taken into account by using a weighted fitting procedure. The fits are acceptable, but 
for high concentrations other fit functions cannot be excluded. 

The concentration dependence of the spin stiffness is given in Fig. (||) and the estimated prefactor C^(S) in Fig. (Q). 
For the pure system (5 = 0) we find 

p s = 0.1734 ±0.0003J 

Q = 0.393 ± 0.008 . (7) 

The shown errors are those which come only from our jsffiighted fitting procedure and are definitely too small. The 
values are in the range of other authors as discussed iraiil. The prefactor slightly increases with concentration up to 
S = 0.25. The spin stiffness decreases by more than a factor of two in this-regime and saturates afterwards. This 
behavior does not agree-^pth the scaling properties proposed by Manousakisu, who used a smaller parameter regime. 

Recently Yanagisawau-3 has predicted making use of the non-linearcr-model with impurities that the exponential 
dependence on the temperature Eq. (Q) changes to a linear dependence 1/T at a very small impurity concentration 
and that £ is a temperature independent constant for S larger than some critical concentration 5 C where the system 
is not ordered. The present data do not quantitatively agree with the work although we cannot rule out a transition 
to the linear dependence. 
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IV. CONCENTRATION DEPENDENCE OF THE LONG RANGE ORDER 



In this chapter the data for the long range order in the ground state are studied. The order parameter of the system 
is the sublattice magnetization 



ieA jeB 



where A and B denote the two different sublattices. In Fig. the sublattice magnetization is plotted against the 
temperature for different hole concentrations and for the system size L — 8. Fig. (||b) shows the data for L = 12 and 
Fig. (||c) for L = 16. For L — 8 there is no big variation for the low temperatures and the extrapolation to T = is 
smooth, which means that the system is near the ground state at the low temperatures. But for L = 12 and L = 16 
the errorbars are enlarged at low temperatures, although we could estimate a general tendency of the concentration 
dependence. In order to obtain these data, up to 9 • 10 5 MCS have been performed for lattices with several Trotter 
numbers. Still the errorbars are rather large. This shows the general difficulty to obtain low temperature data with 
QMC. 

We regard the data at the lowest temperature as the ones in the ground state and extrapolate these data to the 
thermodynamic limit, namely L to infinity. The extrapolation is done using a polynomial in 1/L. This is the same 
extrapolation scheme as hxl, where it is assumed that the spin correlation function decays by a power law c(r) ~ - as 
is indicated by the spin wave theory. In Fig. (|^) the data are plotted against 1/L for some concentrations. The size 
dependences are smoothly extrapolated as well as the pure case. Two different extrapolations are shown, a linear one 
and one with a polynomial up to the quadratic order. The extrapolated value for 6 = 0.375 is less than zero. That 
means the extrapolation scheme is no longer applicable here and the order parameter would already vanish. 

To compare our sublattice magnetization with the spontaneous magnetization rw = (Sf) calculated in the spin 
wave theory, one has to consider the rotational invariancc in spin space and the difference between S and a . In the 
pure case (6 = 0) we find 



m f = y^(NI) =0.315 ±0.04. (9) 

There are many results of this quantity as are seen in review articles00, e. g. Andersonlll found in spin wave 
calculations m? = 0.303 and Reger and Youngo found m) = 0.30 ± 0.02 with QMC. Our result is a little bit higher. 

In Fig. (Q), the linear extrapolated data and the data at T = 0.1 for various sizes are plotted against the concentra- 
tion. The concentration dependences for each size are nearly described by a straight line—This suggests that there is 
only a small cooperative effect of the holes for the global order. The local effect found inB seems to be averaged out. 
The error bars are large in the middle concentration regime. The deviation of data in different samples is also large 
here. This observation implies that the order parameter depends more on the special arrangement of the holes for these 
concentrations. Here we see the extrapolated long rage order decreases with the concentration. It should be noted 
that the concentration dependence has a concave shape in contrast to a convex shape of the concentration dependence 
of the number of the percolating cluster which is equivalent to the long range order in the classical model. This 
suggests that quantum fluctuations have a stronger effect at higher concentrations. In the theory of the non-linear-cr 
model, it is suggested that there is a critical value for the quantum fluctuations. If the system has smaller quantum 
fluctuations than the critical value, the quantum fluctuations are irrelevant and the classical picture of ordering is 
realized. For example, the pure case belongs to this case and also cases with low concentrations obviously belong to 
this case. If the concentration becomes large, however, quantum fluctuations might exceed the critical value. Then 
the system shows a kind of a quantum disordered state where the long range order disappears although the lattice 
itself percolates. This critical value may correspond to S c in the referenceEH and then his conclusion qualitatively 
agree with the present results. 

Although the fluctuation is rather large for definite conclusion, we would estimate the lowest concentration where 
the order parameter becomes zero in Fig. (0). We interpolate the data points and find roughly 

5 C = 0.345 ± 0.015 . (10) 

The error is only interpolated from the errorbars of the extrapolated data points and might be too small. In order to 
obtain a better value of the critical concentration, it is necessary to simulate larger systems at lower temperatures near 
our critical concentration and to average over more different hole configurations. But at present this would exceed 
the ability of standard computers. However rthe above estimation implies that the critical concentration is lower than 
the classical percolation threshold S c = 0.41B. This point should be studied more in the future. 
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There are only very few studies on the doped case. Pimentel and OrbachcS find in spin wave calculations to the 
t—J model S c = 0.32, which they interpret as an upper bound for the sublattice magnetization. But the order in their 
model is destroyed by the mobility of the impurities. 



V. SUMMARY AND DISCUSSION 



We have extended the usual QMC algorithm to treat inhomogeneities in a lattice. Our new flip types enabled us 
to simulate diluted Heisenberg antiferromagnets on the square lattice. Results on the correlation length and the long 
range order have been presented. 

From the spin correlation function, temperature and concentration dependences of the correlation lengths are 
investigated. The correlation lengths show an exponential dependence on 1/T over a wide concentration regime. The 
spin stiffness decreases more than 50% when the impurity concentration increases from zero to 0.25. It saturates 
around S — 0.25 and it seems to stay finite even at concentrations where the sublattice magnetization vanishes. But 
data for high concentrations have a rather large deviation and we do not exclude that the temperature dependence 
of the correlation length changes from the exponential form to a linear form at high concentrations. 

The sublattice magnetization at T = is estimated by extrapolation of data with finite L to infinity. For the pure 
case we get the same results of other authors. The sublattice magnetization of doped systems shows a similar size 
dependence as the pure case. We find a monotonic decrease of the order parameter with hole concentration. From the 
global shape of the concentration dependence, we conclude that the quantum fluctuations have a stronger effect at 
the high concentration region. Although the error bar is large, the critical concentration where the system becomes 
disordered is estimated as S c = 0.345. This value is smaller than the classical percolation threshold. It is suggested 
that quantum fluctuation causes a quantum disordered state in the ground state. 

It is difficult to compare directly the present data with the experimental data for doped LaiCuO^ because of 
the following two reasons. First, in the two-dimensional Heisenberg model the Mermin- Wagner theoremEil forbids 
breaking of a continuous symmetry in two-dimensional models at finite temperatures, while in the materials there is a 
small intcrplanc coupling which creates a non zero Neel-temperature. Second, in the .materials, doping creates mobile 
holes which cause different effects from the static holes. But recently Cheong et. al.@ investigated La^CuxsMgsO^ 
and La2Cui-$ZnsOi, where doping leads to non-magnetic static impurities, which is closer to the present model 
although they are still itinerant systems. They find that the Neel-temperature goes to zero for a concentration 
5 C = 0.22 and that the critical concentration is clearly different from materials which can be described by a model 
with S = 5/2 or by an Ising-model. So there is also a hint for a reduction of the classical percolation threshold from 
experiments. 

We hope that the concentration dependence of the long range order which is presented in this paper is found 
experimentally in some real materials and furthermore that the problem, whether the 8 C and the percolation threshold 
is different or not, is solved in future. This seems difficult for numerical methods at the present situation. We also 
hope for the appearance of very high performance special computers which may solve the problem. 

The authors would like to thank Professor Mikeska for his continuous encouragement. The numerical simulations 
for the present study was very extensive. We spent more than 1500 CPU-hours on the FUJITSU S400/40 (5 Gflops 
peak performance) of the Regionales Rechenzentrum Niedersachsen and several thousand hours on IBM RS6000 
workstations. The present study has been partially supported by a Grant-in- Aid for scientific Research on Priority 
Area from Ministry of Education, Science and Culture of Japan. 
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FIG. 1. The concentration dependence of the correlation function at T = 0.1. 
FIG. 2. The temperature dependence of the correlation length. 
FIG. 3. The concentration dependence of the correlation parameter ps- 
FIG. 4. The concentration dependence of the correlation parameter C^. 
FIG. 5. Temperature-dependence of (AT 2 ) for (a) L = 8 (b) L = 12 and (c) L = 16. 
FIG. 6. Size dependences of (iV 2 ) 
FIG. 7. Concentration-dependence of (iV 2 ) for different system sizes and extrapolated values. 
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